On this page

Spectral Clustering

Spectral Clustering

Overview

The current notebook will demonstrate a simplified machine learning approach to observe the change in a lake water’s extent across time. In order to identify the water, we can use spectral clustering to classify each grid cell into a category based on the similarity of the combined set of pixels across wavelength-bands in our image stacks.

Our example approach uses a version of spectral clustering from dask_ml that is a scalable equivalent of what is available in scikit-learn. We will begin this approach with a single image stack and then conduct a direct comparison on the results from different time points.

This workflow uses data from Microsoft Planetary Computer but it can be adapted to work with any data ingestion approach from this cookbook.

Prerequisites

Concepts

Importance

Notes

Data Ingestion - Planetary Computer

Necessary

scikit-learn

Helpful

Spectral clustering

dask_ml

Helpful

Spectral clustering at scale

  • Time to learn: 20 minutes.

Imports

# Data
import numpy as np
import odc.stac
import pandas as pd
import planetary_computer
import pystac_client
import xarray as xr
from dask.distributed import Client
from pystac.extensions.eo import EOExtension as eo

# Analysis
from dask_ml.cluster import SpectralClustering

# Viz
import hvplot.xarray

Loading Data

Let’s start by loading some Landsat data. These steps are covered in the Data Ingestion - Planetary Computer prerequisite.

Search the catalog

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

bbox = [-118.89, 38.54, -118.57, 38.84]  # Region over a lake in Nevada, USA
datetime = "2017-06-01/2017-09-30"  # Summer months of 2017
collection = "landsat-c2-l2"
platform = "landsat-8"
cloudy_less_than = 1  # percent

search = catalog.search(
    collections=["landsat-c2-l2"],
    bbox=bbox,
    datetime=datetime,
    query={"eo:cloud_cover": {"lt": cloudy_less_than}, "platform": {"in": [platform]}},
)
items = search.get_all_items()
print(f"Returned {len(items)} Items:")
[[i, item.id] for i, item in enumerate(items)]

/usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/pystac_client/item_search.py:841: FutureWarning: get_all_items() is deprecated, use item_collection() instead.
  warnings.warn(

Returned 3 Items:

[[0, 'LC08_L2SP_042033_20170718_02_T1'],
 [1, 'LC08_L2SP_042033_20170702_02_T1'],
 [2, 'LC08_L2SP_042033_20170616_02_T1']]

Load a dataset

item = items[1]  # select one of the results

assets = []
for _, asset in item.assets.items():
    try:
        assets.append(asset.extra_fields["eo:bands"][0])
    except:
        pass

cols_ordered = [
    "common_name",
    "description",
    "name",
    "center_wavelength",
    "full_width_half_max",
]
bands = pd.DataFrame.from_dict(assets)[cols_ordered]
bands
common_name description name center_wavelength full_width_half_max
0 red Visible red OLI_B4 0.65 0.04
1 blue Visible blue OLI_B2 0.48 0.06
2 green Visible green OLI_B3 0.56 0.06
3 nir08 Near infrared OLI_B5 0.87 0.03
4 lwir11 Long-wave infrared TIRS_B10 10.90 0.59
5 swir16 Short-wave infrared OLI_B6 1.61 0.09
6 swir22 Short-wave infrared OLI_B7 2.20 0.19
7 coastal Coastal/Aerosol OLI_B1 0.44 0.02 
ds_2017 = odc.stac.stac_load(
    [item],
    bands=bands.common_name.values,
    bbox=bbox,
    chunks={},  # <-- use Dask
).isel(time=0)

Retain CRS Attribute

epsg = item.properties["proj:epsg"]
ds_2017.attrs["crs"] = f"epsg:{epsg}"

da_2017 = ds_2017.to_array(dim="band")
da_2017

<xarray.DataArray (band: 8, y: 1128, x: 949)>
dask.array<stack, shape=(8, 1128, 949), dtype=uint16, chunksize=(1, 1128, 949), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
  * x            (x) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.637e+05
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
  * band         (band) object 'red' 'blue' 'green' ... 'swir22' 'coastal'
Attributes:
    crs:      epsg:32611

Reshaping Data

The shape of our data is currently n_bands, n_y, n_x. In order for dask-ml / scikit-learn to consume our data, we’ll need to reshape our image stacks into n_samples, n_features, where n_features is the number of wavelength-bands and n_samples is the total number of pixels in each wavelength-band image. Essentially, we’ll be creating a vector of pixels out of each image, where each pixel has multiple features (bands), but the ordering of the pixels is no longer relevant to the computation.

By using xarray methods to flatten the data, we can keep track of the coordinate labels ‘x’ and ‘y’ along the way. This means that we have the ability to reshape back to our original array at any time with no information loss!

flattened_xda = da_2017.stack(z=("x", "y"))  # flatten each band
flattened_t_xda = flattened_xda.transpose("z", "band")
flattened_t_xda

<xarray.DataArray (z: 1070472, band: 8)>
dask.array<transpose, shape=(1070472, 8), dtype=uint16, chunksize=(1070472, 1), chunktype=numpy.ndarray>
Coordinates:
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
  * band         (band) object 'red' 'blue' 'green' ... 'swir22' 'coastal'
  * z            (z) object MultiIndex
  * x            (z) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.637e+05
  * y            (z) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
Attributes:
    crs:      epsg:32611

Standardize Data

Now that we have the data in the correct shape, let’s standardize (or rescale) the values of the data. We do this to get all the flattened image vectors onto a common scale while preserving the differences in the ranges of values. Again, we’ll demonstrate doing this first in NumPy and then xarray.

with xr.set_options(keep_attrs=True):
    rescaled_xda = (flattened_t_xda - flattened_t_xda.mean()) / flattened_t_xda.std()
rescaled_xda

<xarray.DataArray (z: 1070472, band: 8)>
dask.array<truediv, shape=(1070472, 8), dtype=float64, chunksize=(1070472, 1), chunktype=numpy.ndarray>
Coordinates:
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
  * band         (band) object 'red' 'blue' 'green' ... 'swir22' 'coastal'
  * z            (z) object MultiIndex
  * x            (z) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.637e+05
  * y            (z) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
Attributes:
    crs:      epsg:32611

Info

Above, we are using a context manager “with xr.set_options(keep_attrs=True):” to retain the array’s attributes through the operations. That is, we want any metadata like ‘crs’ to stay with our result so we can use ‘geo=True’ in our plotting.

As rescaled_xda is still a Dask object, if we wanted to actually run the rescaling at this point (provided that all the data can fit into memory), we would use rescaled_xda.compute().

ML pipeline

Now that our data is in the proper shape and value range, we are ready to conduct spectral clustering. Here we will use a version of spectral clustering from dask_ml that is a scalable equivalent to operations from Scikit-learn that cluster pixels based on similarity (across all wavelength-bands, which makes it spectral clustering by spectra!)

client = Client(processes=False)
client

Client

Client-cdd263ca-2268-11ee-8ce2-000d3a3b0d78

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: http://10.1.0.50:8787/status

Cluster Info

Now we will compute and persist the rescaled data to feed into the ML pipeline. Notice that our X matrix below has the shape: n_samples, n_features as discussed earlier.

X = client.persist(rescaled_xda)
X.shape

(1070472, 8)

First we will set up the model with the number of clusters, and other options.

clf = SpectralClustering(
    n_clusters=4,
    random_state=0,
    gamma=None,
    kmeans_params={"init_max_iter": 5},
    persist_embedding=True,
)

This next step is the slow part. We’ll fit the model to our matrix X. Depending on your setup, it could take seconds to minutes to run depending on the size of our data.

%time clf.fit(X)

/usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/dask/base.py:1437: UserWarning: Running on a single-machine scheduler when a distributed client is active might lead to unexpected results.
  warnings.warn(

CPU times: user 60 s, sys: 11.5 s, total: 1min 11s
Wall time: 1min 14s

SpectralClustering(gamma=None, kmeans_params={'init_max_iter': 5}, n_clusters=4,
                   persist_embedding=True, random_state=0)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Let’s check the shape of the result:

labels = clf.assign_labels_.labels_.compute()
labels.shape

(1070472,)

labels

array([3, 3, 1, ..., 1, 1, 1], dtype=int32)

The result is a single vector of cluster labels.

Un-flattening

Once the computation is done, we can use the coordinates of our input array to restack our output array back into an image. Again, one of the main benefits of using xarray for this stacking and unstacking is that it keeps track of the coordinate information for us.

Since the original array is n_samples by n_features (90000, 6) and the cluster label output is (90000,), we just need the coordinates from one of the original features in the shape of n_samples. We can just copy the coordinates from the first input feature and populate is with our output data:

template = flattened_t_xda[:, 0]
output_array = template.copy(data=labels)
output_array

<xarray.DataArray (z: 1070472)>
array([3, 3, 1, ..., 1, 1, 1], dtype=int32)
Coordinates:
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
    band         <U3 'red'
  * z            (z) object MultiIndex
  * x            (z) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.637e+05
  * y            (z) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
Attributes:
    crs:      epsg:32611

With this new output array with coordinates copied from the input array, we can unstack back to the original x and y image dimensions by just using .unstack().

unstacked_2017 = output_array.unstack()
unstacked_2017

<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast

<xarray.DataArray (x: 949, y: 1128)>
array([[3, 3, 1, ..., 3, 3, 3],
       [3, 3, 3, ..., 3, 3, 3],
       [3, 1, 3, ..., 3, 3, 3],
       ...,
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1]], dtype=int32)
Coordinates:
  * x            (x) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.637e+05
  * y            (y) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
    band         <U3 'red'
Attributes:
    crs:      epsg:32611

Finally, we can visualize the results! By hovering over the resulting imge, we can see that the lake water has been clustered with a certain label or ‘value’.

raw_plot_2017 = da_2017.sel(band="red").hvplot.image(
    x="x", y="y", geo=True, xlabel="lon", ylabel="lat", datashade=True, cmap="greys", title="Raw Image 2017",
)

result_plot_2017 = unstacked_2017.hvplot(
    x="x", y="y", cmap="Set3", geo=True, xlabel="lon", ylabel="lat", colorbar=False,  title="Spectral Clustering 2017",
)

raw_plot_2017 + result_plot_2017

/usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/operation/__init__.py:14: HoloviewsDeprecationWarning: 'ResamplingOperation' is deprecated and will be removed in version 1.17, use 'ResampleOperation2D' instead.
  from holoviews.operation.datashader import (

Spectral Clustering for 1988

We have conducted the spectral clustering for 2017 and now we want to compare this result to the lake in 1988. Let’s load data from 1988 and run the same analysis as above.

We will use the same catalog, but we will search it for a different point in time and different Landsat mission

Load the data

bbox = [-118.89, 38.54, -118.57, 38.84]  # Region over a lake in Nevada, USA
datetime = "1988-06-01/1988-09-30"  # Summer months of 1988
collection = "landsat-c2-l2"
platform = "landsat-5"  # Searching through an earlier landsat mission
cloudy_less_than = 1  # percent

search = catalog.search(
    collections=["landsat-c2-l2"],
    bbox=bbox,
    datetime=datetime,
    query={"eo:cloud_cover": {"lt": cloudy_less_than}, "platform": {"in": [platform]}},
)

items = search.get_all_items()
item = items[1]  # select one of the results

/usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/pystac_client/item_search.py:841: FutureWarning: get_all_items() is deprecated, use item_collection() instead.
  warnings.warn(

Notice that Landsat 5 data from 1988 has slightly different spectra than Landsat 8 from 2017. Details like this are important to keep in mind when performing analyses that directly compare across missions.

assets = []
for _, asset in item.assets.items():
    try:
        assets.append(asset.extra_fields["eo:bands"][0])
    except:
        pass

cols_ordered = [
    "common_name",
    "description",
    "name",
    "center_wavelength",
    "full_width_half_max",
]
bands = pd.DataFrame.from_dict(assets)[cols_ordered]
bands
common_name description name center_wavelength full_width_half_max
0 red Visible red TM_B3 0.66 0.06
1 blue Visible blue TM_B1 0.49 0.07
2 lwir Long-wave infrared TM_B6 11.45 2.10
3 green Visible green TM_B2 0.56 0.08
4 nir08 Near infrared TM_B4 0.83 0.14
5 swir16 Short-wave infrared TM_B5 1.65 0.20
6 swir22 Short-wave infrared TM_B7 2.22 0.27 
ds_1988 = odc.stac.stac_load(
    [item],
    bands=bands.common_name.values,
    bbox=bbox,
    chunks={},  # <-- use Dask
).isel(time=0)

epsg = item.properties["proj:epsg"]
ds_1988.attrs["crs"] = f"epsg:{epsg}"

da_1988 = ds_1988.to_array(dim="band")
da_1988

<xarray.DataArray (band: 7, y: 1128, x: 949)>
dask.array<stack, shape=(7, 1128, 949), dtype=uint16, chunksize=(1, 1128, 949), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
  * x            (x) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.637e+05
    spatial_ref  int32 32611
    time         datetime64[ns] 1988-07-02T18:03:59.010013
  * band         (band) object 'red' 'blue' 'lwir' ... 'nir08' 'swir16' 'swir22'
Attributes:
    crs:      epsg:32611

Reshape and Standardize

flattened_xda = da_1988.stack(z=("x", "y"))
flattened_t_xda = flattened_xda.transpose("z", "band")
with xr.set_options(keep_attrs=True):
    rescaled_xda = (flattened_t_xda - flattened_t_xda.mean()) / flattened_t_xda.std()
rescaled_xda

<xarray.DataArray (z: 1070472, band: 7)>
dask.array<truediv, shape=(1070472, 7), dtype=float64, chunksize=(1070472, 1), chunktype=numpy.ndarray>
Coordinates:
    spatial_ref  int32 32611
    time         datetime64[ns] 1988-07-02T18:03:59.010013
  * band         (band) object 'red' 'blue' 'lwir' ... 'nir08' 'swir16' 'swir22'
  * z            (z) object MultiIndex
  * x            (z) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.637e+05
  * y            (z) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
Attributes:
    crs:      epsg:32611

Spectral Clustering

X = client.persist(rescaled_xda)
clf = SpectralClustering(
    n_clusters=4,
    random_state=0,
    gamma=None,
    kmeans_params={"init_max_iter": 5},
    persist_embedding=True,
)

%time clf.fit(X)

/usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/dask/base.py:1437: UserWarning: Running on a single-machine scheduler when a distributed client is active might lead to unexpected results.
  warnings.warn(

CPU times: user 37.6 s, sys: 5.24 s, total: 42.8 s
Wall time: 52.5 s

SpectralClustering(gamma=None, kmeans_params={'init_max_iter': 5}, n_clusters=4,
                   persist_embedding=True, random_state=0)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
labels = clf.assign_labels_.labels_.compute()
labels.shape

(1070472,)

labels

array([3, 3, 3, ..., 3, 3, 3], dtype=int32)

Unstack and Visualize

template = flattened_t_xda[:, 0]
output_array = template.copy(data=labels)
unstacked_1988 = output_array.unstack()

<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast

unstacked_1988

<xarray.DataArray (x: 949, y: 1128)>
array([[3, 3, 3, ..., 1, 1, 1],
       [3, 3, 3, ..., 1, 1, 1],
       [3, 3, 3, ..., 1, 1, 1],
       ...,
       [1, 1, 1, ..., 3, 3, 3],
       [1, 1, 1, ..., 3, 3, 3],
       [1, 1, 1, ..., 3, 3, 3]], dtype=int32)
Coordinates:
  * x            (x) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.637e+05
  * y            (y) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
    spatial_ref  int32 32611
    time         datetime64[ns] 1988-07-02T18:03:59.010013
    band         <U3 'red'
Attributes:
    crs:      epsg:32611
raw_plot_1988 = da_1988.sel(band="red").hvplot.image(
    x="x", y="y", geo=True, xlabel="lon", ylabel="lat", datashade=True, cmap="greys", title="Raw 1988"
)

result_plot_1988 = unstacked_1988.hvplot(
    x="x", y="y", cmap="Set3", geo=True, xlabel="lon", ylabel="lat", colorbar=False, title="Spectral Clustering 1988",
)

raw_plot_1988 + result_plot_1988

Spectral Clustering Over Time

Our hypothesis is that the lake’s area is receding over time and so we want to visualize the potential change. Let’s first visually compare the plot of the clustering results from the different time points.

result_plot_1988 + result_plot_2017

By hovering over the lake in the center of each image, we can see that the water was labeled with a cluster label value in both images. Let’s programmatically grab the cluster label at the center for each image.

def get_center_value(arr):
    center_y = arr.shape[0] // 2
    center_x = arr.shape[1] // 2
    center_value = arr[center_y, center_x]
    return int(center_value.values)

water_cluster_1988_label = get_center_value(result_plot_1988.data.value)
water_cluster_2017_label = get_center_value(result_plot_2017.data.value)

Now, we want to align the cluster label for water. Let’s try to set anything that is water to 1 and otherwise 0.

water_1988 = xr.where(unstacked_1988 == water_cluster_1988_label, 1, 0)
water_2017 = xr.where(unstacked_2017 == water_cluster_2017_label, 1, 0)

water_1988.hvplot(geo=True)

water_1988_plot = water_1988.hvplot(
    x="x", y="y", cmap="greys", geo=True, colorbar=False, title="1988 Water"
)

water_2017_plot = water_2017.hvplot(
    x="x", y="y", cmap="greys", geo=True, colorbar=False, title="2017 Water"
)

water_1988_plot + water_2017_plot

Now we can take the difference of these water label arrays to see exactly where the water levels has changed.

with xr.set_options(keep_attrs=True):
    water_diff = water_1988 - water_2017

Red pixels (array value ‘1’) of our image below are where water was lost from 1988 to 2017.

water_diff.hvplot(
    x="x", y="y", cmap='coolwarm', geo=True, xlabel="long", ylabel="lat", colorbar=False, title="Water Change 1988-2017",
)

We did it! We are observing the change in the lake shoreline over time using a simple spectral clustering approach.

Let’s finish things off by adding some geo tiles as a background. To only display the colored pixels overlaid on geo tiles, we could either set the array’s background value (‘0’) to ‘Not a Number’ (NaN), or we could just inform hvPlot that we want the background valued pixels to be transparent with .redim.nodata(value=0).

water_diff.hvplot(
        x="x", y="y", width=400, height=400, cmap='coolwarm', geo=True, xlabel="lon", ylabel="lat", alpha=1, colorbar=False, title="Water Loss from 1988 to 2017", tiles="ESRI",
).redim.nodata(value=0)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/IPython/core/formatters.py:974, in MimeBundleFormatter.__call__(self, obj, include, exclude)
    971     method = get_real_method(obj, self.print_method)
    973     if method is not None:
--> 974         return method(include=include, exclude=exclude)
    975     return None
    976 else:

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/dimension.py:1290, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1283 def _repr_mimebundle_(self, include=None, exclude=None):
   1284     """
   1285     Resolves the class hierarchy for the class rendering the
   1286     object using any display hooks registered on Store.display
   1287     hooks.  The output of all registered display_hooks is then
   1288     combined and returned.
   1289     """
-> 1290     return Store.render(self)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/options.py:1425, in Store.render(cls, obj)
   1423 data, metadata = {}, {}
   1424 for hook in hooks:
-> 1425     ret = hook(obj)
   1426     if ret is None:
   1427         continue

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:280, in pprint_display(obj)
    278 if not ip.display_formatter.formatters['text/plain'].pprint:
    279     return None
--> 280 return display(obj, raw_output=True)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:248, in display(obj, raw_output, **kwargs)
    246 elif isinstance(obj, (CompositeOverlay, ViewableElement)):
    247     with option_state(obj):
--> 248         output = element_display(obj)
    249 elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    250     with option_state(obj):

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:142, in display_hook.<locals>.wrapped(element)
    140 try:
    141     max_frames = OutputSettings.options['max_frames']
--> 142     mimebundle = fn(element, max_frames=max_frames)
    143     if mimebundle is None:
    144         return {}, {}

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:188, in element_display(element, max_frames)
    185 if type(element) not in Store.registry[backend]:
    186     return None
--> 188 return render(element)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:69, in render(obj, **kwargs)
     66 if renderer.fig == 'pdf':
     67     renderer = renderer.instance(fig='png')
---> 69 return renderer.components(obj, **kwargs)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:399, in Renderer.components(self, obj, fmt, comm, **kwargs)
    396 embed = (not (dynamic or streams or self.widget_mode == 'live') or config.embed)
    398 if embed or config.comms == 'default':
--> 399     return self._render_panel(plot, embed, comm)
    400 return self._render_ipywidget(plot)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:406, in Renderer._render_panel(self, plot, embed, comm)
    404 doc = Document()
    405 with config.set(embed=embed):
--> 406     model = plot.layout._render_model(doc, comm)
    407 if embed:
    408     return render_model(model, comm)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:509, in Renderable._render_model(self, doc, comm)
    507 if comm is None:
    508     comm = state._comm_manager.get_server_comm()
--> 509 model = self.get_root(doc, comm)
    511 if config.embed:
    512     embed_state(self, model, doc,
    513                 json=config.embed_json,
    514                 json_prefix=config.embed_json_prefix,
    515                 save_path=config.embed_save_path,
    516                 load_path=config.embed_load_path,
    517                 progress=False)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:560, in Renderable.get_root(self, doc, comm, preprocess)
    543 """
    544 Returns the root model and applies pre-processing hooks
    545 
   (...)
    557 Returns the bokeh model corresponding to this panel object
    558 """
    559 doc = init_doc(doc)
--> 560 root = self._get_model(doc, comm=comm)
    561 if preprocess:
    562     self._preprocess(root)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:146, in Panel._get_model(self, doc, root, parent, comm)
    144 if root is None:
    145     root = model
--> 146 objects = self._get_objects(model, [], doc, root, comm)
    147 props = dict(self._init_params(), objects=objects)
    148 model.update(**self._process_param_change(props))

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:131, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    129 else:
    130     try:
--> 131         child = pane._get_model(doc, root, model, comm)
    132     except RerenderError:
    133         return self._get_objects(model, current_objects[:i], doc, root, comm)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:367, in HoloViews._get_model(self, doc, root, parent, comm)
    365     plot = self.object
    366 else:
--> 367     plot = self._render(doc, comm, root)
    369 plot.pane = self
    370 backend = plot.renderer.backend

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:442, in HoloViews._render(self, doc, comm, root)
    439     if comm:
    440         kwargs['comm'] = comm
--> 442 return renderer.get_plot(self.object, **kwargs)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/renderer.py:70, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     63 @bothmethod
     64 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     65     """
     66     Given a HoloViews Viewable return a corresponding plot instance.
     67     Allows supplying a document attach the plot to, useful when
     68     combining the bokeh model with another plot.
     69     """
---> 70     plot = super().get_plot(obj, doc, renderer, **kwargs)
     71     if plot.document is None:
     72         plot.document = Document() if self_or_cls.notebook_context else curdoc()

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:241, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    238     defaults = [kd.default for kd in plot.dimensions]
    239     init_key = tuple(v if d is None else d for v, d in
    240                      zip(plot.keys[0], defaults))
--> 241     plot.update(init_key)
    242 else:
    243     plot = obj

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/plot.py:936, in DimensionedPlot.update(self, key)
    934 def update(self, key):
    935     if len(self) == 1 and ((key == 0) or (key == self.keys[0])) and not self.drawn:
--> 936         return self.initialize_plot()
    937     item = self.__getitem__(key)
    938     self.traverse(lambda x: setattr(x, '_updated', True))

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/bokeh/plot.py:107, in GeoPlot.initialize_plot(self, ranges, plot, plots, source)
    105 def initialize_plot(self, ranges=None, plot=None, plots=None, source=None):
    106     opts = {} if isinstance(self, HvOverlayPlot) else {'source': source}
--> 107     fig = super(GeoPlot, self).initialize_plot(ranges, plot, plots, **opts)
    108     if self.geographic and self.show_bounds and not self.overlaid:
    109         from . import GeoShapePlot

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:2493, in OverlayPlot.initialize_plot(self, ranges, plot, plots)
   2491 self.tabs = self.tabs or any(isinstance(sp, TablePlot) for sp in self.subplots.values())
   2492 if plot is None and not self.tabs and not self.batched:
-> 2493     plot = self._init_plot(key, element, ranges=ranges, plots=plots)
   2494     self._init_axes(plot)
   2495 self.handles['plot'] = plot

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:487, in ElementPlot._init_plot(self, key, element, plots, ranges)
    480 """
    481 Initializes Bokeh figure to draw Element into and sets basic
    482 figure and axis attributes including axes types, labels,
    483 titles and plot height and width.
    484 """
    485 subplots = list(self.subplots.values()) if self.subplots else []
--> 487 axis_types, labels, plot_ranges = self._axes_props(plots, subplots, element, ranges)
    488 xlabel, ylabel, _ = labels
    489 x_axis_type, y_axis_type = axis_types

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:398, in ElementPlot._axes_props(self, plots, subplots, element, ranges)
    396 # Get the Element that determines the range and get_extents
    397 range_el = el if self.batched and not isinstance(self, OverlayPlot) else element
--> 398 l, b, r, t = self.get_extents(range_el, ranges)
    399 if self.invert_axes:
    400     l, b, r, t = b, l, t, r

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:67, in ProjectionPlot.get_extents(self, element, ranges, range_type)
     65     (x0, x1), (y0, y1) = proj.x_limits, proj.y_limits
     66     return (x0, y0, x1, y1)
---> 67 extents = super(ProjectionPlot, self).get_extents(element, ranges, range_type)
     68 if not getattr(element, 'crs', None) or not self.geographic:
     69     return extents

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/plot.py:1799, in GenericOverlayPlot.get_extents(self, overlay, ranges, range_type)
   1798 def get_extents(self, overlay, ranges, range_type='combined'):
-> 1799     subplot_extents = self._get_subplot_extents(overlay, ranges, range_type)
   1800     zrange = isinstance(self.projection, str) and self.projection == '3d'
   1801     extents = {k: util.max_extents(rs, zrange) for k, rs in subplot_extents.items()}

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/plot.py:1794, in GenericOverlayPlot._get_subplot_extents(self, overlay, ranges, range_type)
   1792         sp_ranges = util.match_spec(layer, ranges) if ranges else {}
   1793     for rt in extents:
-> 1794         extent = subplot.get_extents(layer, sp_ranges, range_type=rt)
   1795         extents[rt].append(extent)
   1796 return extents

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:73, in ProjectionPlot.get_extents(self, element, ranges, range_type)
     71     extents = None
     72 else:
---> 73     extents = project_extents(extents, element.crs, proj)
     74 return (np.NaN,)*4 if not extents else extents

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/util.py:103, in project_extents(extents, src_proj, dest_proj, tol)
    101     geom_in_src_proj = geom_clipped_to_dest_proj
    102 try:
--> 103     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    104 except ValueError:
    105     src_name =type(src_proj).__name__

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:808, in Projection.project_geometry(self, geometry, src_crs)
    806 if not method_name:
    807     raise ValueError(f'Unsupported geometry type {geom_type!r}')
--> 808 return getattr(self, method_name)(geometry, src_crs)

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:944, in Projection._project_polygon(self, polygon, src_crs)
    942     is_ccw = True
    943 else:
--> 944     is_ccw = polygon.exterior.is_ccw
    945 # Project the polygon exterior/interior rings.
    946 # Each source ring will result in either a ring, or one or more
    947 # lines.
    948 rings = []

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/geometry/polygon.py:99, in LinearRing.is_ccw(self)
     96 @property
     97 def is_ccw(self):
     98     """True is the ring is oriented counter clock-wise"""
---> 99     return bool(self.impl['is_ccw'](self))

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:14, in is_ccw_impl.<locals>.is_ccw_op(ring)
     13 def is_ccw_op(ring):
---> 14     return signed_area(ring) >= 0.0

File /usr/share/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:7, in signed_area(ring)
      3 """Return the signed area enclosed by a ring in linear time using the
      4 algorithm at: https://web.archive.org/web/20080209143651/http://cgafaq.info:80/wiki/Polygon_Area
      5 """
      6 xs, ys = ring.coords.xy
----> 7 xs.append(xs[1])
      8 ys.append(ys[1])
      9 return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(ring.coords)))/2.0

IndexError: array index out of range

:Overlay
   .Tiles.I :Tiles   [x,y]
   .Image.I :Image   [x,y]   (value)

Summary

Starting from raw Landsat data, we have used a simple spectral clustering approach to observe the change in a lake water’s extent across time.

What’s next?

Adapt this notebook for your own use case or select another workflow example notebook.

Resources and References

  • Authored by Demetris Roumis circa Jan, 2023

previous

Preprocessing - Regrid